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EXTENDING THE MODEL OF KH 15D: ESTIMATING THE EFFECTS OF FORWARD SCATTERING AND 

CURVATURE OF THE OCCULTING RING EDGE 

Devin W. Silvia 1 ' 2 and Eric Agol 1 

To appear in the Astrophysical Journal 

ABSTRACT 

The periodic eclipses of the pre-main-sequence binary, KH 15D, have been explained by a circumbi- 
nary dust ring inclined to the orbital plane, which causes occultations of the stars as they pass behind 
the ring edge. We compute the extinction and forward scattering of light by the edge of the dust ring 
to explain (1) the gradual slope directly preceding total eclipse, (2) the gradual decline at the end of 
ingress, and (3) the slight rise in flux at mid-eclipse. The size of the forward scattering halo indicates 
that the dust grains have a radius of a ~6 (D/3 AU) /an, where D is the distance of the edge of the 
ring from the system barycenter. This dust size estimate agrees well with estimates of the dust grain 
size from polarimetry, adding to the evidence that the ring lies at several AU. Finally, the ratio of the 
fluxes inside and outside eclipse independently indicates that the ring lies at a few astronomical units. 

Subject headings: stars: individual (KH 15D) — stars: pre-main-sequence — circumstellar matter 

ible star is covered. This model provides an excellent 
fit to the characteristics of the KH 15D light curve and 
radial velocities, as shown by Winn et al. (2004) and 
Chiang & Murray-Clay (2004). Chiang k Murray-Clay 
(2004) also suggested that the ring must be somewhat 
warped in order to maintain the nodal precession. This 
picture proved correct when KH 15D was confirmed to be 
an eccentric spectroscopic binary with a two-year, high- 
resolution, multi-site spectroscopic study (Johnson et al. 

2004) . Archival data showed evidence for the second star 
which is completely hidden in recent data (Johnson et al. 

2005) . 

Winn et al. (2006) made a further refinement of the KH 
I5D model by adding a faint blue halo around each star 
within the binary providing an explanation for the grad- 
ual slope of the light curve directly before and just after 
total eclipse, as well as the slight increase in flux dur- 
ing mid-eclipse. In their model, they solved for the one 
dimensional halo brightness with 4 free parameters and 
an arbitrary functional form that resulted in a strongly 
asymmetric shape around both stars due to the asym- 
metric shape of ingress and egress, resulting in a puzzling 
physical picture: why should identical asymmetric halos 
surround both stars? Winn et al. (2006) also included a 
gradual change in the angle of the edge of the ring pro- 
jected onto the sky to account for longer-term variations 
the the lightcurve. 

In this paper, we focus on fitting the photometric data 
for KH 15D spanning 1995 to 2004 and the 18 published 
radial velocity data points. We compute the effects of 
(1) a gradual change in opacity at the ring edge (prior 
models treat the ring edge as being sharp) and (2) for- 
ward scattering by dust at the ring edge, which together 
naturally explain the asymmetric shape of ingress and 
egress (Winn et al. 2006 acknowledge that forward scat- 
tering may be a more sensible explanation for their halo) . 
We do not include the rotation of the edge of the ring 
(Winn et al. 2006), but instead propose that the edge 
of the ring is curved as a means of explaining additional 
variation in the eclipse durations that is not explained 
by the precession of the ring. We show that these mod- 
ifications provide a much better fit to the photometric 



1. INTRODUCTION 

The object KH 15D (V582 Mon) is a binary weak- 
lined T Tauri star in the cluster NGC 2264. As a class, 
T Tauri stars are thought to be young stars that are still 
accreting gas from the remains of their parent molecular 
clouds. This period of stellar evolution lasts for a few 
million years and is characterized by optical variability 
and chromospheric emission lines (see Bertout 1989 for 
an extended discussion). 

Kearns & Herbst (1998) reported a periodic flux vari- 
ation in KH 15D: it drops by 3.5 magnitudes every 48 
days. Such large flux variations (>2 mag in V) are seen 
in "Type III" irregular variables in young clusters (Joy 
1945, Herbst et al. 1994), but the remarkable property 
of KH 15D is that is also periodic. This depth of eclipse 
can not be explained solely as a binary due to the large 
magnitude and duration, so the earliest theories of KH 
15D postulated that circumstellar material eclipsed the 
star. 

Since Kearns & Herbst (1998) first reported the un- 
usual properties of KH 15D, theoretical explanations for 
its properties have included an edge-on circumstellar disk 
(Hamilton et al. 2001, Agol et al. 2003, Winn et al. 2003), 
an asymmetric surrounding envelope (Grinin & Tam- 
bovtseva 2002), and an orbiting vortex of solid particles 
(Barge & Viton 2003). However, it was soon realized 
that KH 15D was likely an eccentric binary system that 
is being occulted by the edge of a circumbinary dust ring 
(Winn et al. 2004, Chiang & Murray-Clay 2004), and 
currently only one star is visible during part of its orbit 
causing the large change in magnitude. In addition, the 
advance of the screen as a function of time due to the 
nodal precession of the ring, which is inclined relative to 
the binary, explains the lengthening of the duration of 
the eclipses as a larger portion of the orbit of the vis- 

1 Department of Astronomy, University of Washington, Box 
351580, Seattle, WA 98105; dwsl29@astro.washington.edu, 
agol@astro.washington.edu 

2 Current address: Department of Astrophysical and Planetary 
Sciences, University of Colorado, UCB 391, Boulder, CO 80309; 
devin.silvia@colorado.edu 



2 



Silvia & Agol 



data of KH 15D and we speculate on what additional in- 
sight this provides us about young stellar evolution and 
protoplanetary processes. In section [5] we summarize the 
data which we use to fit our model. In section [3] we dis- 
cuss the elements included in our model. In section |4] we 
discuss the best-fit model parameters and errors, and in 
section [5] we discuss the implications of these results and 
possible future directions. 

2. DATA 
2.1. Photometry 

We used the 6694 I-band photometric data points from 
Hamilton et al. (2005) taken between 1995 to 2004 with 
a dozen different telescopes using CCDs. We ignored 
older, sparser data which is not precise nor dense enough 
to show the effects of forward scattering. We converted 
the I-band magnitudes to fluxes using /„/ fo = lo' -0 47 ' 
where f ~ 2600 Jy is the flux zero-point. We enlarged 
the error bars for all data points in which the star was 
completely eclipsed or completely uneclipsed due to ob- 
served fluctuations that are larger than the reported er- 
rors. Although these fluctuations may be partly due 
to starspots (Hamilton et al. 2005), for the purposes of 
model fitting we treat these fluctuations as an additional 
source of systematic error and replaced the errors bars 
that accompanied the data with the standard deviation 
of the scatter of the data in or out of eclipse, respectively. 

In addition to the orbital elements of the binary and 
the parameters describing the ring (described below), we 
introduce three primary flux parameters of our model: 
{fin, A, out, f2,out}, where /i„ is the non- variable flux of 
the system, /i iOU t is the flux received when only Star 1 
is visible, and fi,out is the flux received when only Star 
2 is visible. 

2.2. Radial Velocities 

We fit our model to the published radial velocity data 
from Hamilton et al. (2003) and Johnson et al. (2004). 
We did not discard data points taken near eclipse as 
these data points did not appear to show the Rosssiter- 
McLaughlin effect (Rossiter 1924, McLaughlin 1924), so 
we use all 18 values for the data fitting process. We 
inflate the error bars for these radial velocity values to 
account for a systematic radial velocity offset due to light 
from the binary that is scattered off the back and sides 
of the circumbinary ring, an effect discussed by Herbst 
et al. (2008). If we assume that the parameter fi n is 
due to the large angle scattering by the dust ring, which 
will have a range of Doppler shifts between —K and K, 
where K is the velocity semi-amplitude, then the system- 
atic offset in the measured radial velocity values should 
be of order V% (fin/ fi.out) K ■ This offset was added in 
quadrature to the uncertainties published by Winn et al. 
(2006). Our decision to modify the error bars was to 
reduce the discrepancy between the orbital parameters 
found from the radial velocity data versus those found 
from the light curve data, as pointed out by Winn et 
al. (2006). We implement weighting of the contribution 
to x 2 from photometric and radial velocity in the same 
manner as Winn et al. (2006), described in more detail 
in§H 

3. THE KH 15D MODEL 




X (AU) 

Fig. 1. — TOP - The nearly edge-on orbits of the two stars in KH 
15D as would be seen from the observer's perspective. The solid 
curved line represents the edge of the occulting screen when the 
curvature parameter is added to the model while the dashed line 
represents the edge of the occulting screen with zero curvature. 
/3 is the angle between the X-axis and the edge of the screen as 
described in the text. BOTTOM - The orbits of KH 15D as they 
would be seen from above. In this representation, the occulting 
screen would not be visible. 

3.1. Orbit Model 

The binary orbit was computed by solving Kepler's 
equation, and converting the orbital elements to Carte- 
sian coordinates and velocities, following the procedures 
described in Murray and Dermott (1999). The param- 
eters describing the binary orbit are {P, e, i, to, T p , 7}, 
where P is the period, e is the eccentricity, i is the incli- 
nation of the system, w is the argument of pericenter, T p 
is the time of pericenter passage, and 7 is the heliocen- 
tric radial velocity of the center of mass. The geometry 
of the orbit is shown in Figures [1] and [H where the origin 
is the barycenter of the binary, the Z-direction is along 
the line of sight, and we fix the longitude of the ascend- 
ing node at 180° so that the orbits cross the sky plane on 
the X-axis, with the visible star (which we call Star 1) 
crossing the sky plane on the positive X-axis as it moves 
towards positive Z. We initially fit just the radial veloc- 
ity data, confirming that our derived model parameters 
agreed with Winn et al. (2004). 

3.2. Ring Orientation, Occultation, and Curvature 

We define the angle between the ring edge, projected 
onto the X — Y plane, and the X-axis as (3. We allow the 
edge of the ring to move in the Y-direction with velocity 
v r i n g and define the time at which the edge passes the 
center of mass of the system as t ring . Figure Q] shows 
the orbits as they would be seen from the observer's per- 
spective and a snapshot of the position of the edge of the 
ring. To help visualize KH 15D, the approximate three- 
dimensional geometry is of the system is shown in Figure 

m 

We assume the stars are limb-darkened with a 
quadratic limb-darkening law, 

7(r) = l- 7l (l- e )- 72 (l-e) 2 , (1) 

where e is the cosine of the angle between the normal to 
the stellar surface and the line of sight to the observer. 
We set 71 = 0.4478 and 72 = 0.2091, appropriate for a 
star with T eff w 4300K (Agol et al. 2004) in the I-band 
according to the tables of Claret (2001). 

Initially we treat the ring as a "knife edge" which 
changes abruptly from transparent to opaque, and thus 
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Fig. 2. — A three dimensional representation of the geometry of 
the KH 15D system. Although the orbits are defined by the best fit 
model, the distance between the barycenter of the system and the 
circumbinary ring is not fully known and so the plotted distance 
is just one possible value. The shape and extent of the ring is not 
derived from our model, but is only shown to indicate one possible 
geometry. 

any part of the binary system seen above the edge is un- 
obscured, while below the edge is completely occulted. 
For this sharp edge model, the flux during ingress and 
egress is given by 



F(z) = ! + (/) " 



11 



1 



1 - 7i ^"72 + ^72-z' 



1 - 7i - -72 ) cos" 



--7r(7i + 2 72 )(2 - 3z + z 3 ) 
6 



(2) 



where z is the distance between the center of either star 
and the edge of the occulting ring in the units of the 
stellar radius and (I) = tt(1 — ^71 — ^72)- This equation 
is only valid for \z\ < 1 and takes on F(z) — for z < — 1 
and F(z) = 1 for z > 1. Although the knife-edge model 
provides a fairly accurate set of initial parameters, the 
quantitative agreement with the data is poor, giving a 
X 2 of ^12991 for 6697 degrees of freedom, where 6697 = 
6694 photometric points + 18 velocity points - 15 free 
parameters. Most of the discrepancy in this fit occurs 
at the beginning of ingress and end of egress where the 
model has a much steeper slope than the data as shown 
in Figure [3] 

The next step in our calculation is to add "fuzziness" 
to the edge of the ring by choosing the optical depth to 
vary as a power-law as a function of the sky-plane dis- 
tance from the ring edge (it has infinite optical depth 
below the edge). This power-law behavior seems phys- 
ically plausible if we assume that there is a roll-off in 
density at the edge of the ring or due to the vertical 
structure of the ring. We compute the transmitted flux 
during ingress and egress by convolving the knife-edge 
with the power-law optical depth, 



F* = 



max(0,z-\-l) 



x J max(0,z — 1) 



dxe-t*°W (—) 1+a F(z-x) 
i) v x ) 



(3) 

where xq is the scale length of the fuzzy edge (projected 
on the sky) in units of stellar radius, R s tar (we take both 




z (in units of R st „) 

Fig. 3. — The occultation curve for the currently visible star as a 
function of the distance between the center of the star and the point 
on the edge of the occulting screen where the optical depth becomes 
infinite. The dots represent the data, but that data has been binned 
to allow for a clearer view of the agreement between the model and 
the observed values. The solid line represents the full model while 
each of the other lines represents the shape of the curve due to 
different parts or versions of the model. The shape of the forward 
scattering curve is directly dependent on model parameters; the 
height of the curve is determined by K as while the width of the 
curve is determined by w. 

stars to have the same radius); a is the power-law expo- 
nent of optical depth variation; and x is perpendicular 
to the ring edge. In this model, we initially used a fixed 
value of a = 2 but after establishing a reasonable fit to 
the data we allowed a to vary freely. We report w, which 
is xq transformed to physical units, w = xoR star (in me- 
ters). This addition reduces x 2 to 12890 for 6695 degrees 
of freedom, improving the fit by A% 2 = 101. 

We next add curvature of the ring edge to the model 
which improves the fit to mid-eclipse data. We parame- 



terize this as y r 



[IX 



ring 



where y ring and x ring are 



the y and x coordinates perpendicular and parallel to 
the edge of the ring at the origin when it crosses the 
barycenter of the system. We find that by adding cur- 
vature the rate at which the central re-brightcnings fade 
during 1995-1998 is slower than without it since the hid- 
den star's orbit remains visible for a slightly longer period 
of time, agreeing better with the lightcurve data. Again 
we see an incremental improvement in the model as \ 2 
goes to 9582 with 6694 degrees of freedom, producing 
A\ 2 = 3308. We neglected precession of the ring as the 
model of Winn et al. (2006) indicates the edge projected 
on the sky plane would have rotated only ~4.5 degrees 
over the 9 year data set. 

3.3. Forward Scattering 

Forward scattering by dust is the most important new 
addition to our model: the same dust at the edge of the 
ring causing extinction will also diffract light, leading 
to an apparent halo around the star. The same effect 
can be seen on nights when the moon passes behind a 
thin layer of clouds. The water droplets in the clouds 
diffract light from the moon producing a halo of forward- 
scattered light about the moon. For both the moon and 
KH 15D the halo is present when the optical depth is near 
unity — at very high optical depths multiple scatterings 
will cause the radiation to eventually be absorbed, while 
at very low optical depths there is almost no scattering. 
As the star passes behind the edge of the ring, this for- 
ward scattering halo softens the shape of the ingress and 
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egress, as observed in KH 15D. 

The angular distribution of scattered light for spherical 
grains can be modeled approximately as an Airy disk. 
Since the wings of the Airy disk are much weaker than the 
central peak, we approximate the Airy disk by a Gaussian 
angular distribution to allow for faster computation of 
the scattered flux, 



where dap / dQ is the differential scattering cross section 
into a solid angle d£l. A Gaussian with a standard devi- 
ation of ag = 0.43A/(27ra), where a is the radius of the 
scattering particle, gives a good approximation to the 
peak of an Airy disk. There may be additional large- 
angle scattering by the ring; however, if the size of the 
ring is much larger than the orbital size of the binary, 
then this radiation should have a much weaker time de- 
pendence than the forward-scattered light, so we treat 
this as a constant flux contribution during the eclipse, 
namely fi n mentioned in § 2.1. 

As the optical depth becomes large, multiple forward 
scatterings can occur. We make the assumption that the 
scattering region (due to the ring surrounding the two 
stars) can be treated in the thin-screen approximation; 
that is, (1) the scattering angle is small enough that the 
photon can be approximated as going straight as it passes 
through the dust ring, cos (9 scat ) ~ 1, and thus the total 
optical depth for each forward-scattered photon is sim- 
ply the line-of-sight optical depth to the star; and (2) 
the dust is at approximately a constant distance from 
the source. Since absorption and large-angle scattering 
can occur in addition to forward scattering, we add an 
additional parameter, K as , the ratio of the sum of the 
absorption and large-angle scattering opacity to the for- 
ward scattering opacity. 

With these approximations, we can compute the 
specific intensity of scattered radiation, I(x, y) = 

2jv=T In(x, y) at a position (x, y) relative to the source 
position (xo,yo)j both positions projected on the sky 
plane and relative to the ring edge (which is rotated rel- 
ative to the X — Y plane). For radiation that has under- 
gone N scatterings, the specific intensity of the scattered 
radiation, In(x, y), is given by a Gaussian, with standard 
deviation \^Nag, times the probability of scattering N 
times, 

I N (x,y) = 

Z7> n-N(~ ,.\ (*-so) 2 + («-to) 2 
■gO (X, V) r -r(x,y)(l+ Kas ) to^WTl 

2irNa 2 e Nl 

(5) 

where Fq is the flux of the star at the distance of the dust 
scattering screen, D is the star-dust separation, and r is 
the forward-scattering optical depth through the slab at 
position (x,y). We treat the stars as point sources for 
the purposes of computing the forward-scattered light. 
When K as ~ 1, truncating the sum at N max = 5 scat- 
terings gives a converged distribution of scattered light. 
We report <jd = D x ag in meters rather than ag. 

Figure [4] shows the distribution of scattered light as a 
function of distance of the center of the star from the 
edge (t = oo ) of the ring. The best-fit parameters (see 
below) for the KH 15D system were used in creating this 



Fig. 4. — Star and forward scattering halo as a function of the 
height above the edge of the ring. Each frame is 12 stellar radii 
on a side, and y represents the height of the center of the star 
above the edge of the ring (where r = oo) in units of the stellar 
radius. The grey (or color) scale represents the logarithmic I-band 
intensity where white is the surface brightness of the center of the 
star while black is < 10 -4 of this brightness. Limb-darkening of 
the star is included, but the star appears nearly uniform due to 
the logarithmic intensity scale. An mpeg version of this figure is 
available on-line. 

figure. When the star is far above the ring edge (y = 6.0 
frame), the forward scattering is weak as the optical 
depth is small within the angular size of the scattering 
halo. As the star approaches the edge, forward scatter- 
ing increases as the optical depth increases (y — 3.5, 1.0 
frames), but the stellar light starts to decrease due to 
extinction by the dust. As the star becomes eclipsed, 
the halo can still be seen above the edge, but it is much 
weaker due to extinction of the scattered halo (y = —1.5 
frame). 

For computing the sky-plane separation of the star and 
ring edge, we take into account the curvature of the edge 
of the ring since the distance traveled by the star and the 
edge of the ring is comparable to the radius of curvature. 
However, to simplify the scattering calculation we treat 
the edge as being straight (since the radius of curvature 
is much larger than the edge scale-length) and the op- 
tical depth varies only perpendicular to the ring edge, 
allowing analytic integration of the scattered specific in- 
tensity parallel to the edge (along x), while integration 
perpendicular to the edge is carried out numerically. 

After adding forward scattering, the model light curve 
matches the more gradual slope on either end of the 
eclipses as well as the small rises in flux during mid- 
eclipse. The x 2 reduces to 7162 for 6692 degrees of free- 
dom, an overall improvement of A% 2 = 5829 as compared 
to the knife-edge model with no forward-scattering. In 
particular, we notice that the combination of forward 
scattering with a curved ring edge produces mid-eclipse 
"bumps" in the model that are more pronounced and 
widened, which agrees well with the data. Furthermore, 
we see the noticeable effects of the hidden star decrease as 
the heights of the mid-eclipse bumps diminish over time 
(Figure [3]) ■ This supports the idea that as time passes 
an increasing portion of the orbit of the binary system is 
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obscured by the ring of material. 

As a comparison between our model and the "halo" 
model of Winn et al. (2006), we take the flux data as 
a function of the separation between Star 1 and the 
edge of the ring (as in Figure [3J and run a x 2 mini- 
mization using the halo model on the data points for 
which Star 2 contributes an insignificant flux. Using 
the seven parameters necessary for modeling the occupa- 
tion curve, {R s tar, fin, fi,out, w, a D ,K as , a} for our model 

versus {R star-, fin, fl,out, e lj e 2j £,1, £2} fbx Winn et al. S 

(2006) model, we find that forward scattering produces 
X 2 = 5849.4 and the halo model produces x 2 = 5903.7 
for 5484 degrees of freedom. The difference in \ 2 be- 
tween the two models is A% 2 = 54, indicating that the 
forward scattering model is favored over the halo model. 

4. RESULTS 

The complete model has 18 free parameters, consist- 
ing of the 6 orbital elements, {P, e, i, lo, T p , 7}, 4 pa- 
rameters defining the orientation, shape, and motion of 
the occulting screen, {v r i ng , t r i ng , ft, /^}, and 8 parame- 
ters that modify the structure of the lightcurve, {R s tar, 
fin, fi.out, h,out, w, gd, K as , a}. In addition, we as- 
sume the masses of the stars are Mi = 0.6 ±0.1 M Q and 
M2 = 0.72 ± 0.1 M©, respectively. This is based on the 
theoretical pre-main-sequence evolutionary tracks and 
the mass-luminosity relation indicating that the mass of 
the visible star should be 0.6 ±0.1 M and the mass 
ratio, A/2/M1, should be 1.2 ± 0.1 (Winn et al. 2006). 

For the optimization of these parameters we choose to 
minimize with respect to a \ 2 defined in a similar way 
to that of Winn et al. (2006), such that 

N f / f f \2 N v , _ s2 

2 X / JmodA Jdata,i \ , \ \ / ^modA ^dataA \ 

x =51 <■» ) +A Sl ) • 

(6) 

where Nf = 6694 and N v — 18. The presence of A allows 
us to apply a certain amount of weight to the relatively 
few radial velocity measurements so that equally good 
fits are found to both fluxes and radial velocities. To 

determine what value of A to use for our model, we at- 

2 2 

tempted to find a value that would produce « jj?-, 

where %/ an d xl are the separate non-reduced x 2 val- 
ues for the fluxes and radial velocities respectively. This 
procedure gives A = 9 for our final model. Using that 
value of lambda we find that our best-fit parameter set 
produces x 2 = 6992.6 and x 2 = 18-8. A list of the best- 
fit parameters and errors (described next) is presented 
in Table [TJ In addition, Figures [5] and [5] show the model 
and corresponding data for both the lightcurve and the 
radial velocities. 

To estimate the uncertainties in each of our 18 param- 
eters we followed the procedure of Winn et al. (2006). 
We repeatedly fit our model to artificial data sets with 
random noise added to our best-fit model. To model the 
noise, we took the residuals of the best-fit model, nor- 
malized them by the model flux, randomized them, and 
then added them back to the best-fit model (again scal- 
ing to the flux) to create an artificial data set. Since the 
number of radial velocities points is so small, we instead 
added Gaussian noise with a standard deviation equal 
to the error bars of the radial velocity data (Winn et al. 



TABLE 1 
Model Parameters for KH 15D 



r aramctcr 


hjstimated Value 


U nits 


Mi 


0.6 ± 0.1 


Mq 


M 2 


0.72 ± 0.1 


Mq 


P 


48.359 ± 0.0012 


days 


e 


0.51 ± 0.008 


— 


i 


89.31 ± 0.38 


degrees 


OJ 


5.1 ± 1.3 


degrees 


T P 


2352.72 ± 0.13 


J.D.-2450000 


7 


18.6 ± 1.5 


km s -1 


luring 


58.46 ± 5.95 


m s _1 


tring 


-1289.64 ± 340.23 


J.D.-2450000 


Rstar 


1.2 ± 0.23 




Jin 


(5.075 ±0.123) x 10~ 8 


fo 


fl,out 


(1.59 ±0.03) x 10" 6 


fo 


f2,out 


(3.28 ±0.28) x 10~ 6 


fo 


W 


(1.30 ±0.67) x 10 9 


111 


&D 


(4.35 ± 0.34) x 10 9 


111 




1.91 ± 0.18 




P 


1.18 ± 0.09 


radians 


/-< 


2.00 ± 0.40 


ATJ- 1 


a 


1.53 ± 0.59 





2006). For each artificial data set we also allowed Mi and 
M2 to vary with Gaussian distributions that had means 
of 0.6 and 0.72 respectively and standard deviations of 
0.1 so that the resulting error bars would have the ±0.1 
Mq constraint placed on them as mentioned above. To 
determine the error bars on the remaining parameters, 
which arc quoted in Table [T] we took the standard devi- 
ations of 2000 parameter sets that were produced from 
fits to 2000 random realizations of the noisy lightcurve. 
We found that subsequent iterations did not appear to 
produce significant changes in the errors. 

5. SUMMARY AND DISCUSSION 

We find that our best-fit model is in good agreement 
with the data with a reduced x 2 of order unity. The 
model accurately reproduces the fall-off in the emergence 
of Star 2 during the 1995-1998 time frame as well as the 
gradual deepening of the total eclipses from 1998-2004 
(Figure [S]). Furthermore, we find one of the most satis- 
fying results to be the success of forward scattering as 
the source for the mid-eclipse bumps in the lightcurve 
once Star 2 no longer passes beyond the edge of the ring, 
as well as explaining the gradual fall off and rises in the 
curve as Star 1 enters and exits the eclipses (Figure [3]). 
Our model also predicts that around the beginning of 
2008 the photosphere of Star 1 will no longer move be- 
yond the edge of the occulting screen - the same time 
frame as predicted by Winn et al. (2006). After that 
point, the maximum light we will receive from KH 15D 
will be due to forward scattering as Star 1 nears the edge 
of the ring, which will last roughly 5-7 years. 

Several of the predicted parameters are in relatively 
good agreement with the values we expect. First, 7 is 
in agreement with both the value published by Winn et 
al. (2006) as well as the median heliocentric radial veloc- 
ity of the cluster NGC 2264, 20 ± 3 km s" 1 (Soderblom 
et al. 1999). Second, we allowed the radius of Star 1, 
Rstar, to be a free parameter determined by the shape 
of the lightcurve, finding a best-fit of 1.2±O.23i?0 which 
is consistent with the value used in Winn et al. (2006), 
who held the parameter fixed at 1.3 Rq based on the 
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Fig. 5. — The lightcurve of KH 15D as a function of phase. Each plot shows one year worth of data overplottcd with the best fit model. 
The most notable features are the drop off in mid-eclipse rebrightenings from 1995-1998, the gradual deepening of the eclipses from start 
to finish, and the small rises in flux during mid-eclipse from 1998-2004, which are produced by forward scattering in our model as the ring 
edge moves across the orbit of the binary. 

of this parameter. 

We now compare our best-fit model with the model 
in Winn et al. (2006). Our best-fit period, eccentricity, 
and inclination are relatively close to those presented by 
Winn et al. (2006), agreeing to within 0.1%, 10%, and 
5%, respectively. However, formally there is quantita- 
tive disagreement (> la) between Winn et al.'s (2006) 
values and our own which most likely results from phys- 
ical differences between their model and ours, as well as 
fitting different data sets with different assumed errors. 
Our best-fit flux for Star 2 is nearly twice that of Star 1, 
somewhat larger than the ratio of 1.36 found by Winn et 
al. (2006). 

From the forward scattering width in our model, ag, 
we can estimate the average size of the dust grains 
which dominate the forward-scattering opacity. Using 
the parameter a& and our Gaussian approximation ag = 
0.43A/(27ra) with A = 8140 A, we find that the radius 
a of the dust grains is ~6(D/3 AU) yum, where D is the 
distance between the edge of the ring and the stars. This 
agrees extremely well with the ~6-8 [im size of the dust 
grains estimated by Agol et al. (2004) based on the weak 
variation of the observed polarization with wavelength 
(smaller/larger grains cause a stronger/ weaker variation 
of polarization with wavelength than observed). A dis- 
tance of 3 AU was estimated by Chiang & Murray-Clay 
(2004) and Winn et al. (2004) to explain a rate of preces- 
sion that matches the observed velocity of the ring edge 
across the binary, indicating internal consistency between 
the forward-scattering, polarization, and precession con- 
straints on a and D. 



Fig. 6. — The radial velocity profile of KH 15D as a function of 
phase. The solid line represents the model while the dots represent 
the 18 measured radial velocities and their associated errors bars 
as defined in the text. 



mass-radius relation as predicted by pre-main-sequence 
evolutionary tracks. Finally, we find that our value 
for the eccentricity, 0.51, is consistent with the pseudo- 
synchronization upper limit of 0.66 as discussed by Winn 
et al. (2006). 

Our derived value for the angle between the edge of the 
ring and the horizontal axis of the system, (3 = 68° ± 5°, 
was not what we expected. Although our error bars seem 
to indicate that the value is relatively well-constrained, 
since the system is nearly edge on, changes in (3 will not 
produce significant changes in the model lightcurve, so 
we are not confident in the best-fit value or uncertainty 
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Beyond our confirmation that forward scattering is the 
correct physical process to produce the lightcurve fea- 
tures, we can also use the flux parameters to estimate 
the apparent area of the ring. One possible interpreta- 
tion of the residual flux in eclipse is that it is due to 
large-angle scattering off of the ring. This implies that 
the ratio of fi n to the sum of /i, ut and /2,o«t is approx- 
imately the ratio of scattered flux from the ring to the 
total flux of the stars, and thus gives an indication of 
the solid angle covered by the ring. From our param- 
eters, -, — % w O.Of. If we then make a simplify- 

Jl,outiJ2,out 

ing assumption that the ring of occulting material has 
roughly a shape of a warped ring with mean radius R 
and maximum warp angle ip (the difference between the 
inner and outer edges), then the solid angle covered by 
the warped ring is on average ~ 2Trip. Let us then say 
that perhaps we only see the back half of the warped 
ring, that only half of the light that hits the ring is re- 
flected back towards us, and that of order half of the 
back of the ring is obscured by the front of the ring. 
This would imply that f in w ^i/)(fi tOUt + h.out)- Since 
we observe a value of fyj w 0.01, this indicates 

Jl, out +J2, out 

that ip w 0.08. We also know that the warped ring at 
the very least covers the orbit of the stars, and since our 
fit implies that the edge is at a large angle to steller or- 
bits (which are nearly edge-one), which implies Rip > 0.3 
AU and thus R > 0.3/0.08 ~ 4 AU. This argument con- 
tains several highly uncertain geometry-dependent and 
scattering-dependent factors, but is another consistency 
check, independent of the above arguments, that the ring 
lies at a few astronomical units. Further modeling of the 
ring is definitely warranted. 

Additional tests of our model can also be carried out 
to confirm its accuracy. One such test would be to see 
if a scattering halo is present at infrared wavelengths. 
At longer wavelengths, there will be competing effects 
between a larger forward scattering angle and a smaller 
absorption opacity. We are currently applying our model 
to Spitzer data which will be published once data analysis 
is complete. Another test would be to compare the radial 
velocities of the forward-scattered stellar light and the 
large-angle scattered light. The forward-scattered light 
should have a radial velocity equal to that of the eclipsing 
star, while the large-angle scattered light should have a 
wide range of radial velocities as the light is scattered off 
of different parts of the ring (Herbst et al. 2008). 

Our particular choice for the parameterization of the 
optical depth of the ring-edge as a power-law should be 
taken to be fairly generic since for the purposes of mod- 
eling the scattering halo we are only concerned with the 
opacity near r <~ 0.1 — 3 and over this small range of r 
a power-law functional form should provide a good fit to 
almost any smooth variation. Predictions for the varia- 
tion in surface density near a gap caused by a planet (if 
the ring is shepherded by a planet) obey a more complex 
functional form, but depending on how the surface den- 



Agol, E., Barth A. J., Wolf, S., & Charbonneau, D. 2004, ApJ, 600, 
781 

Bertout, C. 1989, ARA&A, 27, 351 

Barge, P. & Viton, M. 2003, ApJL, 593, L117 

Barsunova, O. Yu., Grinin, V. P., Sergeev, S. G. 2005, Astrophysics, 
48, 1 



sity maps to opacity in the I-band, only the small region 
near r ~ 1 can be approximated by a power law. 

The parameter \x is essentially an estimate of the sky- 
plane radius of curvature, \i ~ 1/(2R). Our best-fit value 
of n = 2 AU -1 implies a sky-plane radius of curvature of 
0.25 AU, which is a rather small value if the ring edge is 
located at ~ 3 AU. This might imply strong warping of 
the ring, or it may simply be that this is not the correct 
parameterization of the properties of the ring edge. 

Rice et al. (2006) produced simulations of a ring with 
an imbedded planet around a T-Tauri star. They found 
that the surface density of their rings varies by an order 
of magnitude over a distance that is ~ 0.1 of the semi- 
major axis. In our model the same variation occurs over 
<~ 5w ~ 6 x 10 9 m, about 1/8 of the value predicted by 
the Rice et al. (2006) model, which agrees reasonably well 
if the ring is nearly edge-on causing a smaller sky-plane 
scale-length. This is not at all evidence that a planet is 
present in the system, but a shepherding planet may be 
one explanation for the variation in optical at the edge 
of the ring. The scale length should depend on opacity 
which varies with wavelength, an additional reason to 
make infrared observations. 

To finish our discussion of model parameters, we com- 
ment on n as <~ 2. As the wavelength of light is A = 0.8 
/Km and we have estimated the radius of the dust parti- 
cles is a = 6 /im, then the particle scattering parameter 
is x = 27r(j) « 50. In this limit, the forward scattering 
cross section and reflectance (large-angle scattering plus 
absorption) cross section equal the area of the particle, 
so one would expect n as = 1- Therefore, our value of 
K as rj 2 indicates that there may be grains smaller than 
6 fim mixed in with x rj 1. These smaller grains would 
not have as narrow a forward-scattering peak and would 
thus contribute to the large-angle scattering/absorption 
component — resulting in an increase in the value of K as ■ 
Since large-angle scattering causes polarization and for- 
ward scattering does not, a possible test of the model pro- 
posed in this paper would be to measure the degree of po- 
larization as a function of the ratio of forward-scattered 
to large-angle scattered light, which changes during the 
ingress/egress. Infrared, radial-velocity, and polarimet- 
ric measurements, combined with more detail modeling 
of the dust ring, should help to put tighter constraints on 
the properties of the system and unlock further mysteries 
of KH 15D. 
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